Bounding the mass of the graviton with gravitational waves: Effect of spin precessions 

in massive black hole binaries 
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Observations of gravitational waves from massive binary black hole systems at cosmological dis- 
tances can be used to search for a dependence of the speed of propagation of the waves on wavelength, 
and thereby to bound the mass of a hypothetical graviton. We study the effects of precessions of 
the spins of the black holes and of the orbital angular momentum on the process of parameter 
estimation based on the method of matched filtering of gravitational-wave signals vs. theoretical 
template waveforms. For the proposed space interferometer LISA, we show that precessions, and 
the accompanying modulations of the gravitational waveforms, are effective in breaking degeneracies 
among the parameters being estimated, and effectively restore the achievable graviton-mass bounds 
to levels obtainable from binary inspirals without spin. For spinning, precessing binary black hole 
systems of equal masses 10^ Mq at 3 Gpc, the lower bounds on the graviton Compton wavelength 
achievable are of the order of 5 x 10^® km. 

PACS numbers; 



I. INTRODUCTION AND SUMMARY 

The anticipated launch of the Laser Interferometer 
Space Antenna (LISA) in the 2020 time frame will pro- 
vide a promising new tool for doing astrophysics with 
massive binary black hole systems. The inspiral and 
merger of massive black holes (MBH) with masses of the 
order of 10^ - 10^ Mq will be detectable to large distances 
in LISA'S sensitive frequency band between 10^^ and 1 
Hz. The detection of gravitational waves (GW) from 
MBH systems will allow us to infer important astrophys- 
ical and astronomical information, such as the masses 
and spins of the black holes, the location of the system 
on the sky and its distance from the solar system. 

Another important aspect of MBH binaries is the pos- 
sibility of testing general relativity itself. In previous pa- 
pers we have studied the bounds that could be placed on 
alternative theories of gravity such as scalar-tensor theo- 
ries of the Brans-Dicke type, and theories in which grav- 
itational waves propagate with a wavelength-dependent 
speed, as if the "graviton" were massive [l[ [3, S, 13, H, @] . 
Specifically, in [sj we showed that the inclusion of aligned, 
non-precessing spins weakens the bounds obtainable on 
the graviton mass by almost an order of magnitude. This 
is because the parameters that characterize the inpiral- 
ing binary are highly correlated, so that the addition of 
parameters (the spins) into the estimation process effec- 
tively dilutes the available information, leading to weak- 
ened bounds or estimates on most parameters. 

However, Vecchio pointed out that when the effects 
of precession of spins are incorporated into the gravi- 
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tational waveforms, i.e. when the spins are not aligned 
with the orbital angular momentum, the accuracy of pa- 
rameter estimation can be improved. He studied the so 
called "simple precession" case where either one of the 
bodies has zero spin, or the black hole masses are equal, 
and only spin-orbit interactions are included. The mod- 
ulations of the amplitude and phase of the gravitational 
waveform induced by the precession of the spin(s) and 
by the precession of the orbital plane effectively adds in- 
formation to the estimation process, partially decouples 
some of the parameters, and thus leads to restored ac- 
curacy. Lang and Hughes P extended Vecchio's work 
to include arbitrary spins and masses, and also spin-spin 
interactions, and found significant improvements in the 
accuracy of mass measurements as well as sky localiza- 
tion. In addition, they showed that the magnitudes of the 
spins of the binary members, especially for low redshift 
systems at z ~ 1, could be measured with accuracies of 
the order of 10~^. 

In this paper we describe the results of an indepen- 
dent code written by one of us (AS) for analysing binary 
inspiral with precessing spin and for carrying out parame- 
ter estimation based on the method of matched filtering, 
but extended to include the effects of a massive gravi- 
ton. In addition to confirming the central conclusions of 
Lang and Hughes we show that spin precessions sig- 
nificantly improve the bounds that can be placed on the 
mass of the graviton. In parallel work, we have shown 
that including higher signal harmonics in the PN wave- 
form (but without spins) also leads to improved bounds 
on the graviton mass @. 

Our main conclusion, shown in Figs. [1] and [2] is that 
the inclusion of spin precession effects increases the lower 
bound on the graviton Compton wavelength Ag by almost 
an order of magnitude, on average, with respect to the 
one calculated for the same non precessing system. Recall 
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FIG. 1: Distribution of lower bounds on the graviton Comp- 
ton wavelength \g (in units of 10^^ km) for 10* equal-mass 
(10® Mq) black- hole binaries at redshift z = 0.55, or a lumi- 
nosity distance 3 Gpc, randomly located on the sky. Number 
of bins is set to 50. First three histograms (narrow lines; red, 
blue, green in the color version) assume either no spins or 
aligned spins with spin-orbit (SO) and/or spin-spin (SS) cou- 
pling. Final two histograms (thick lines; violet and black in 
the color version) include precessions induced by non-aligned 
spins 



that Xg is related to the mass of the graviton by Xg = 
h/rUgC, where h is Planck's constant and c is the speed 
of light, so that a lower bound on Xg represents an upper 
bound on rrig. 

Indeed, the new bounds, labeled MG+SO+PREC and 
MG-hSO-hSS-hPREC in Fig. [U which incorporate spin- 
orbit (SO) effects only and spin-orbit and spin-spin (SS) 
effects, respectively, along with the effect of a mas- 
sive graviton (MG), are comparable to those inferred 
from an identical system without spin effects at all, la- 
beled MG. This improvement is independent of mass, 
as seen in Fig. [2l which plots median lower bounds 
on Xg for systems without spin (MG), with non pre- 
cessing spins (MG-I-SO-I-SS), and with precessing spins 
(MG+SO+SS+PREC) for various pairs of masses span- 
ning two orders of magnitude in total mass. 

The rest of the paper provides the details of the analy- 
sis behind our main conclusion. In Sec. |TT] we review the 
construction of gravitational waveform templates and the 
orbital dynamics when spin precessions are incorporated. 
In Sec, mil we describe the parameter estimation process 
based on the method of matched filtering, and in Sec. 
ITVlwe present the results. Section IVl presents concluding 
remarks. Throughout the paper we use units in which 
c = G = 1. 



II. GRAVITATIONAL WAVEFORM AND 
ORBITAL DYNAMICS INCLUDING SPIN 
PRECESSIONS 



100 



3 







1 








• 


■ ■ • MG 










■ MG+SO+SS 








♦— 


» MG+SO+SS+PREC 



































/' 










/ 










/ 










« 



































1 







lo' lo' lo' 



Total Mass (M^ ) 

FIG. 2: Median lower bounds on the graviton Compton wave- 
length \g (in units of 10^^ km) for 10* black-hole binaries 
at redshift z = 0.55, or a luminosity distance 3 Gpc, ran- 
domly located on the sky. Systems contain black holes of 
mass (1, 1) X 10^ (1, 10) x 10^ (1, 1) x 10^ (1, 10) X 10® and 
(1, 1) X 10^ Mq, from left to right, respectively. 



In this section we give a brief overview of the assump- 
tions made for the GW signal used for our calculations. 
The waveform emitted by an inspiraling black hole bi- 
nary system can be described accurately by the post- 
Newtonian approximation developed by several groups 
(see for example for a review of the post-Newtonian 
approximation for gravitational wave emission from in- 
spiraling binaries see For our study we made 
the following assumptions, some physically justified and 
some imposed for simplicity, (i) We take into account 
only the inspiral phase of the signal, ignoring the merger 
and ringdown part. The bound on the graviton mass 
will be dominated by information from the inspiral phase 
where the wavelength of the signal varies over many or- 
ders of magnitude, (ii) We assume the restricted second 
post-Newtonian (2PN) approximation, in which the am- 
plitude of the signal is evaluated to the lowest, Newtonian 
order, while the phase is evaluated to 2PN order. (Ref. 
Q goes beyond this approximation, but does not include 
spins.) (iii) We use the stationary phase approximation 
(SPA) for calculating the Fourier transform of the signal, 
(iv) We assume that the orbits are quasi-circular. 

With these assumptions and following Q, we express 
the Fourier transform of the GW signal as 

rT~_-2/3 Ayf5/6 
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where / is the frequency of the wave; M is the "chirp 
mass" of the system given by = rf^^M, with M = 
mi + 1712 and rj — mim2/A/^; 1 — 1,2 labels two 
possible combinations of data from the three arms of 
LISA. The quantity t{f) is the time at which the emitted 
gravitational-wave frequency equals /. The distance of 
the source D^, is given as a function of the redshift by 
the expression, 



Dl{z) 



(1 + ^) 



dz' 



(2.2) 



with the cosmological parameters having the values 
= 0.25, flA = 0.75, Ho = 75kms-i Mpc'^, following 
the latest fits by the WMAP mission [ill • 

The amplitude of the wave, A^[t{f)], is given by the 
expression, 



-4(L • n)2F/ {tf 



1/2 



(2.3) 



where, L and n are unit vectors in the directions of the 
source orbital angular momentum and the line of sight 
to the source, respectively. The LISA antenna pattern 
functions for one pair of arms F^'^ are given by the 
expressions [T2| . 

F+(6's,0s,V's) = ^{1 + cos^ 9s) cos 2(j)s cos 2i;s 
— cos 6*5 sin 2(j)s sin 2ips , 

F^{0s,<l>s,i^s) = ^{l+cos^0s)cos2<j>ssm2il;s 

+ cos 6*5 sin 2(/)s cos 2?/;5 , (2.4) 

where 9$ and (jjs are the spherical angles for the binary's 
line of sight n in a frame in which the three LISA space- 
craft are at rest, and ips is the polarization angle of the 
wave in the same frame given by the expression 



tan -05 



q • z L • z — (L • n) (z • n) 



P z 



n • (L X z) 



(2.5) 



The unit vector z is orthogonal to the plane of the LISA 
satellites, while p and q are axes orthogonal to fi, defined 
as p = n X L/|n x L| and q = p x n; they are the 
principal axes of the wave, i.e. defined such that the two 
polarizations are exactly 90° out of phase. For the second 
"detector" (actually a linear combination of outputs from 
the three LISA arms such that the noise is independent 
of the noise in the two arms that make up detector 1) the 
expressions are, 

F,^es,cj,s,^s) = F,^{es,q^s-^,i^s). (2.6) 

In order to use these expressions for our calculations 
they must be transformed to a coordinate system tied to 



the ecliptic. Taking into account the "cartwheel" motion 
of the LISA array as it orbits the Sun, we use expressions 
in Ref. M, 



cos 9 s 



1 - \/3 

- cos 9s — sin 9s cos[(f>{t) — cjys] 



t 



(fis = ao + 27r— + tan ^ /3 , 

^ ^ \/3cosgg -fsinggcos[(^(t) - (t>s] 
2 sin ^5 sin[^(t) — ^g] ' 

where, 9s and (ps denote the fixed direction to the source, 
and (j){t) = 00 + 2'Kt/T denotes barycentric longitude of 
the detector's center-of-mass as it orbits the Sun, where T 
is one year and cto are arbitary orientation constants 
usually chosen to be zero. The polarization angle "05 is 
written in terms of barycentric angles using Eq. (j2.5[) 
and the expressions 

1 - 

z-ii = -cos 6*5 ^sin6'scos(0(t) - 0s) , 

L-z icos^L - ^sin ^icos (0(t) - 0l) , 

L • n = cos ^lcos 9s 

-hsin^L s\n9s cos (0i - 05) , 

n • (L X z) = i sin 9l sin 9s sin (0l - 0s) 
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cos 9l sin 9s sin(0(t) — 0s) 



-cos^s sin6'L sin(0(i) - 0l)] . (2.8) 

where 9^ and 0l are the polar and azimuthal angles of 
the orbital angular momentum vector L in barycentric 
coordinates. 

In order to determine t{f) in Eq. (|2.ip . we use the 
rate at which the observed frequency changes because 
of the emission of gravitational radiation by the binary 
system and because of the propagation delay induced by 
a massive graviton, as given by the expression [T^j, 
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dt 
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where 



Pa 



n'^DM 



(2.9) 



(2.10) 



describes the contribution of the massive graviton. Its 
effect is to alter the time of arrival of the wavefronts for a 
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given frequency, as a function of the Compton wavelength 
Xg and a distance parameter defined as 



D 



1 + z 



dz' 



Ho Jo (l + z')V["M(l + 2')' + r!A] 



(2.11) 



The parameter /3, describing spin-orbit interactions, is 
given by 



L • S, 



(2.12) 



and the parameter a describing spin-spin interactions is 
given by 



V 



721(L-Si)(L-S2)- 247(Si-S2) 



(2.13) 



where Xi = Si/mf, is the dimensionless spin parameter 
for each body. 



To get the relation between the time elapsed and the 
frequency, one has to integrate Eq. (|2.9|) . When spin pre- 
cessions are taken into account, both the spin-orbit and 
spin-spin coefficients /3 and a are oscillating functions 
of time around an average value; however, as shown in 
[l^ . the amplitude of the oscillations is small so one can, 
without significant loss of acccuracy, assume that they 
are constant for the purpose of the integration. The re- 
sult is 



/ 3058673 5429 
V 1016064 1008 



(2.14) 



In our calculations we use the above expression wherever necessary to express time as a function of frequency, but 
we insert the frequency-dependent values of (3 and a that come out of the numerical integration of the spin precession 
equations, as described below. Although this is a slightly inconsistent procedure, we do not expect it to have a large 
effect, since the spin-orbit and spin-spin terms are high-order PN corrections, and thus are relatively small. 

The phase $/ in Eq. (|2.ip has several terms that describe different effects contributing to the phasing of the 
gravitational wave, in the form 



(2.15) 



The first term ^(/), is the phasing fimction at 2PN order arising from the internal dynamics of the binary system, 
given by the expression |;5|], 



*(/) = 2nft,-^, 



TT 
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128 



-5/3 



743 
336 



4(47r-/3)(7rM/) + 10 



/ 3058673 



V 1016064 1008 



l-^/3.(-A^/r/^ + f 



5429 617 2 \ , ,,.^4/3 
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(2.16) 



r 



where tc and <i>c, are the time and the phase of coales- 
cence respectively. Here, as in the calculation of t(/), we 
hold (3 and a fixed during the required integrations, and 
then insert the time- varying values afterward. 

The term '■pl,o\\t{f)\, often called the "polarization 
phase" , arises from the conversion of the real signal into 
an amplitude p.3p and a phase, and is given by the ex- 
pression. 



<^poi(<) = tan 1 



2(L-n)F/(t) 
[l + (L.n)2]F+(t) 



(2.17) 



The term ¥'D[i(/)], is the "Doppler phase", arising from 
the varying arrival time of the signal as the detector 
moves around the sun, given for both detectors by the 
expression, 

'PD{t) = 27r/i?e sin 0s cos[0(i) - 4>s] , (2.18) 

where i?© = 1 AU. 

Finally, the term Sp^[t{f)], comes from the integrated 
change in the orbital phase, caused by the precession of 
the orbital angular momentum vector that accompanies 
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the spin precessions, and is given by [1 
-/final 2L • n 



Spmf)] 



df 



f 



1 - (L . n) 



(Lxn)-L. (2.19) 



In contrast to the case where the spins are ahgned with 
the orbital angular momentum, precessions of the spins 
and of the orbital plane induce modulations of both the 
amplitude and phase of the gravitational wave on a pre- 
cession time scale. The orbital time scale is given by 



while the precession time scale is given by 

5 ^5/2 

" ^ /iAfi/2 ■ 



precession 



(2.20) 



(2.21) 



The final relevant time scale is that of the inspiral, given 
by 



inspiral 



■T2 



(2.22) 



Since for most of the inspiral, Toibitai < T'prcccssion < 
Tlnspirai, wc are justified to use orbit averaged equations 
for the spin and angular momentum precessions, and to 
allow the total angular momentum J to evolve adiabati- 
cally as a result of gravitatjxjnal radiation damping. 
The relevant equations [isj are 



Si 
L 



where. 



Hi = — 



Oi X Si , 

X S2 , 

J — Si — S2 



2 mi ' 



(2.23a) 
(2.23b) 
(2.23c) 



-^(S2-L)L + is2 




are the orbit-averaged precession vectors, and 



j = -f^(^)^/^L, 
5 r r 



(2.24a) 



(2.24b) 



(2.25) 



is the change in the total angular momentum due to radi- 
ation reaction to lowest PN order and for a quasi-circular 
orbit. The overdot denotes a usual time derivative. 

From Eqs. (|2.23a|) and (j2.23bp . the magnitudes of the 
spin vectors Si do not change, so the dimensionless spin 
parameters Xi constant. When |L| ^ |S|, spin-orbit 
coupling dominates, and the rate of precession of each 



spin is independent of the other spin. In other words, 
when spin-orbit effects dominate, binaries with slowly 
spinning objects produce roughly as many precession cy- 
cles as do binaries with faster spinning objects. The dif- 
ference is that for small |S| the cone describing the pre- 
cession of L is smaller. 

In the general case of arbitrary initial conditions for the 
spin and angular momentum vectors the above system of 
equations cannot be solved analytically, so we must resort 
to numerical integration using routines from [l^. 



III. PARAMETER ESTIMATION 

We carry out the parameter estimation using the 
standard theory of Fisher information matrices and 
the maximum likelihood approximation developed for 
gravitational-wave applications by several authors [3, 
117,,I8,]. 

Given the noise spectrum of the instrument and a sig- 
nal h{t] 0°') characterized by a number of parameters 6'° 
of the source, one can define the inner product between 
two signals hi{t), h2{t) as follows, 



(/ll|/l2) = 2 df 



M{f)h2{f) + hl{f)h,{f) 



= 4Re 



df 



Snif) 
Snif) ' 



(3.1) 



where hi{f) and /i2(/) are the Fourier transforms of the 
respective gravitational waveforms hi(t^ 0°), and the star 
denotes complex conjugate, and Sn{f) is the noise spec- 
tral density of the detector. The signal-to-noise ratio 
(SNR) for a given signal h{t) is then given by, 



p[h] ^ {h\hf/^ 



(3.2) 



evaluated at the estimated values 0°" of the source pa- 
rameters. In our analysis we will include the possibility 
that both detector combinations of LISA will be opera- 
tional. In this case, the Fisher information matrix Tab of 
the source is defined as follows, 



Tafc = 



(3.3) 



where h^, and are the signals in the two LISA arm 
combinations discussed earlier. In the limit of large SNR 
and if the noise is stationary and Gaussian, the probabil- 
ity that the GW signal h{t) is characterized by a given 
set of values of the source parameters 6*° is 



p{0\h) = p("H6/)exp 



1 



(3.4) 



where A6'° — 9°^ — 6°-, is the difference between the esti- 
mated and the true value of the parameters, and p^"-* (9) 
is the prior information. An estimate of the rms error 
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in measuring the source parameter 0° can then be cal- 
culated, in the limit of large SNR, by taking the square 
root of the diagonal elements of the inverse of the Fisher 

matrix, 

A9^^,, = V^, I] = r-i. (3.5) 

Finally, the correlation coefficients between two parame- 
ters 6*" and 9'' are given by 



and for extra-galactic binaries by [24I , 

Sr'^^'if ) = 4.2 X 10-4" Hz-i , 



(3.10) 



Cab 



(3.6) 



It turns out that because LISA is designed to detect 
massive inspirals, will naturally provide the largest lower 



The total noise spectrum to be used [2^ is given by 
Su{f) = min{5r*'-(/)exp(^T„;.i,,„„d7V/d/), 

Sr''{f) + SfU)} + 5r-^^'(/) , (3.11) 

where Tmission is the duration of the mission, which we 
assume to be one year, k = 4.5 is the average number of 
frequency bins that are lost when each galactic binary is 
fitted out, and dN/df = 2 x IQ-^/n/a Hz-i. 

In calculating the integrals for the Fisher matrix, we 



bounds on Xg. This can be seen from the dependence of l^^^^ ^^.^ following egressions for the lower and upper lim- 
theboundonAg on the relevant parameters of the system its of mtegration 91 . 
and the detectors, given by Eq. (4.9) of 0: 



its of integration The initial frequency is given by 



oc 



1/4 



D 



1/2 _yv^ll/12 

il/4 „l/3 ' 



(3.7) 



where Sq is a parameter that establishes the floor of the 
noise spectral density (in Hz"^), /o is a characteristic 
"knee" frequency, or frequency where the noise is a min- 
imum. The quantities /(7) and A are determined from 
the Fisher matrix inversion and are largely independent 
of either 5o or /q, or of the SNR of the signal. In any 
case, the bound is only weakly dependent on these vari- 
ables. The ratio D/{\ + Z)Dl is weakly dependent on 
distance, reflecting the fact that the effect of the massive 
graviton and the estimation errors both grow with dis- 
tance. Finally, the factor f^,^'^ is roughly the same 
for LISA as it is for, say, advanced LIGO, and thus the 
best bound on Ag will come from LISA. 

The noise spectrum of LISA consists of the instrumen- 
tal noise intrinsic to the on-board instrumentation and 
drag-free control, and astrophysical noise due to unre- 
solved astrophysical sources of GW lying in the instru- 
ment's frequency band. The instrumental noise currently 
used in the literature is that of reference [l^ (also found 
online at [20j). In our calculations we use an analytic 
version of the instrumental noise following [s*] , given by. 



^instr^y.) ^ [9.18 X IQ-^V""^ 

+9.18 X 10-^** f] 



I- 1.59 X IQ-^^ 



(3.8) 



where / is in Hz. Technically this model ignores the os- 
cillatory effects in the transfer function of LISA at high 
frequencies where the gravitational wavelength becomes 
comparable to the spacecraft separations, but since the 
relevant systems for bounding the graviton mass are mas- 
sive binary inspirals at the low frequency end, we do not 
expect this simplification to have a large effect. 

The spectral density for the noise from galactic binaries 
is approximated by , 



Sf\f) = 2.1 X lO"'^''^ f-''/^ Hz" 



/initial = niax{/iow, /(Tobs)} 

/(Tobs) = 4.149 X 10 



VlOH/o 



-5/8 



yr 



-3/8 

(3.1?) 



where /(Tobs) comes from the leading term of Eq. (j2.14|) . 
with Tobs = ic — t{f)j and where /low is the lower cutoff 
of the LISA instrument, taken here to be 10^^ Hz. The 
final frequency is given by 



/final — niin{/lSCOi /end} , 



(3.13) 



(3.9) 



where /isco = (6^/^7rM)^^ is the usual frequency for 
the innermost stable circular orbit and /end = 1 Hz is a 
conventional upper cutoff for the LISA noise curve. In 
order to see clearly the effects of spin-precessions on the 
graviton-mass bound, we choose the same observation 
time Tobs = 1 year as in 



IV. RESULTS 

In general a quasi-circular binary black hole inspiral 
in GR is described by a total of 15 parameters; adding 
the parameter for the massive graviton, we have the fol- 
lowing 16 parameters: the two individual masses of the 
system, In(TOi), ln(m2), the luminosity distance to the 
source In(T'i), the two dimensionless spin parameters 
Xi, X2, the time and phase at coalescence tc, <&c, the 
two angles of the binary's sky position (f>s, cos9s, the 
two angles of the initial orientation of the orbital angu- 
lar momentum vector, 0^, cos^l, the four angles of the 
initial orientations of the spins of the two bodies, (j)Si , 
cos^Sj, cos^Sj, and finally /3g, the parameter that 
describes the massive graviton contribution to the phase 
of the waveform. All angles are defined in the frame at- 
tached to solar system barycenter. 

The inclusion of spin precessions makes some of the 
parameters used traditionally for estimation less suitable. 
For example, the spin-orbit and spin-spin parameters /3 
and (T are now time (frequency) dependent, so one must 
either go directly to the values of xi a-nd X2 and the four 
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initial spin orientation angles as parameters, or one must 
use the initial values of /3 and a along with four other 
suitable parameters (such as the initial angles) as the 
appropriate parameters. We choose the former. Instead 
of the chirp mass (M) and the symmetric mass ratio (77), 
we use the individual masses as parameters, because they 
more directly scale the spins. 

Proceeding with the error estimation, we first fix a pair 
of masses in the source rest frame, the phase at coales- 
cence and a redshift or luminosity distance. We then 
randomly select the dimensionless spin parameters, Xi 
within the range [0, 1], and the initial spin, orbital angu- 
lar momentum and source position angles (eight angles) . 
We also select randomly the time of the coalesence tc, 
within the assumed time duration of the mission, which 
corresponds to different orientations of the LISA arms at 
the first reception of the gravitational wave signal. One 
effect of the selection of random tc , is that some signals 
might be partially cut off because they are already in the 
sensitive band when LISA starts observing them. We use 
the routine RAN2 [l^ to produce the random numbers. 

We also set the nominal value /3g = for our calcula- 
tions since we are interested in setting a lower limit for 

The inclusion of spin precessions modulates both the 
amplitude and the phase of the waveform. Since the total 
angular momentum J = Si-|-S2-fLis conserved on a 
precession timescale, the orbital angular momentum vec- 
tor L must precess to cancel out the efi^ects of spin pre- 
cessions. As a consequence, the amplitude given by ex- 
pression ()2.3|) , now changes and modulates the waveform 
accordingly. The phase is also affected mainly through 
the terms that describe the polarization phase (|2.17[) and 
integrated change in orbital phase (|2.19p . Finally, in the 
phasing function ^'(/) the parameters (3 and a are now 
frequency dependent. 

Another thing to note about precession is the follow- 
ing. Since we generate arbitrarily the initial directions 
of the spin and angular momentum vectors we can have 
both kinds of precession, simple and transitional, as de- 
scribed in [1^. "Simple" precession is the (most com- 
mon) case where the angular momentum vector L and 
the total spin vector S precess around the total angu- 
lar momentum vector J, which decreases slowly because 
of gravitational radiation reaction. Simple precession al- 
ways occurs when |L| 3> |S|, which is generally the case 
early in the inspiral. "Transitional" precession occurs 
when L and S are almost antialigned and |L| < |S|. It 
consists of a "tumbling" of the L and S vectors (with 
the sum still tied to J) because of the loss of "gyroscopic 
bearings" of the system. 

Apostolatos et al. found that, in order to get tran- 
sitional precession, the initial angle between the total 
spin S and the angular momentum L must be larger than 
about 164°, so that, as |L| decreases because of radiation 
reaction, the conditions for transitional precession will be 
met during the inspiral phase. We have checked our ini- 
tial values and learned that out of the 10^ sets of initial 
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Bound on A, (10 km) 



FIG. 3: Distribution of lower bounds on the graviton Comp- 
ton wavelength Xg (in units of 10^^ km ) for 10* binaries when 
spin is included without precession. The system is a 10® -I- 
lO'' Mq BBH at z = 0.55 (3 Gpc). Solid (dashed) lines refer 
to one (two) LISA detectors. Number of bins is set to 50. 
Results match Fig. 6 of 5j 



angles for Si, S2 and L only about 80 lead to transitional 
precession, and for these, the Fisher matrix calculations 
were not adversely affected by the complicated preces- 
sions. 

Our calculations start with the numerical integration 
of the spin precession equations (|2.23ap - (|2.23cp in the 
frequency domain, using the random initial values for the 
six parameters of the spins of the two bodies and the two 
components of the orbital angular momentum, to get the 
orientations of L, Si, S2 over the duration of the signal. 
We use Eq. (|2.9p to convert from d/dt to d/df, and use 
Kepler's third law at lowest order, r = M^^"^ /{nf)'^/^ to 
convert from r to / in Eq. (|2.24p . We use a fourth order 
Runge-Kutta constant step size routine RK4 [16]. Once 
this is done, the spin parameters (3 (spin-orbit), a (spin- 
spin) and the integrated phase correction dp^[t{f)] (|2.19[) 
are calculated. Subsequently the signal in the frequency 
domain, Eq. (|2.ip . is calculated on the same grid on which 
the precession equations are solved. 

All the derivatives of the signal with respect to pa- 
rameters needed for the Fisher matrices are calculated 
numerically. Given a determination of h{f) for a given 
set of initial values of the 16 parameters 0", we also cal- 
culate h{f) for nearby values 0° -I- 66'^ and 6^ ~ SO"" for 
each parameter in turn. Then for each 0°, we calculate 
the derivative using the standard central finite difference 
formula, 

(4.1) 
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Bound on ;i (lo'^km) Correlation coefficient 



FIG. 4: Distribution of lower bounds on \g (in units of 10^^ 
km) for 10^ binaries including spin precessions. The system 
is a 10*^ + 10*^ Mq BBH at z ^ 0.55 (3 Gpc). Red curve is 
for tc fixed to one year; blue curve is for random values of tc 
in the one year interval of the LISA mission. Solid (dashed) 
lines refer to one (two) LISA detectors respectively. 



for each value of / on the grid. Since we are using double 
precision accuracy for our variables a natural choice of 
the small shifting parameter 66'^ for the calculation of 
the numerical derivatives would be 56°' ~ 10~^ — 10~*. 
We have chosen S9°' = 10^'* for all the sixteen parameters 
estimated in order to achieve the best possible accuracy. 
Then the necessary integrals are calculated numerically 
on the same grid using the extended Simpson rule for the 
closed interval [/initial, /final]. 

Finally, the inverse of the Fisher matrix is calculated 
using the routine SVDCMP [l^. We have also used 
LU and Gauss- Jordan decomposition as a cross check, 
with identical results. The main advantage of this rou- 
tine is that it allows us to check whether the matrix is 
ill-conditioned for the inversion, by calculating the ra- 
tio of the smallest to the largest eigenvalue of the ma- 
trix. If this ratio is of the order of the machine accu- 
racy (~ 10~^^), then the matrix inversion is not to be 
trusted. However, another simple sanity test is to mul- 
tiply the original with the produced inverse matrix and 
check how far the product is from the identity matrix. 
This can be measured by the maximum deviation of the 
non-diagonal elements of the matrix from zero (~ 10^^^) 
in double precision. We have checked both of these cri- 
teria. In most of the cases the condition number is of 
the order of 10^^'' — 10^^^ and the maximum deviation 
is of the order of 10"^ — 10~^. In the cases where the 
condition number approaches double precision accuracy 
i.e. ~ 10~^^, the maximum deviation is of the order of 
10-4. 




Correlation coefficient 



FIG. 5: Distribution of the correlation coefficients between 
spin parameters Xi s^nd /3g and between individual masses rrn 
and f3g 



We have carried out several tests and diagnostics for 
the validity of our code. In the case of aligned, non- 
precessing spins, we have reproduced the fifth panel of 
Fig. 6 of Berti et al. @] for one year integration time 
and tc fixed; the results, shown in Fig. [3l agree within 
the natural statistics of our Monte Carlo simulations. 
In contrast to [5|, we have been able to quote errors on 
the graviton mass including spin-spin effects because of 
the improved machine precision available that allowed us 
to invert the larger matrices reliably. Also for the non- 
precessing case, and for individual choices of angles, we 
have compared parameter estimation errors and correla- 
tion coefficients with those from a Mathematica code de- 
veloped independently and used by K. G. Arun for other 
calculations; the agreement was excellent. 

For the precessing cases, we have checked our code 
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with respect to the median error results quoted by Lang 
and Hughes [1, f^or the asymmetric mass systems 
of (mi,m2) = (3,1) X IO^Mq and {mi, 1712) = (3,1) x 
IO^Mq at z — 1,3,5 respectively. Modulo the statistics 
of the Monte Carlo simulation, we found good agreement 
for the median errors in masses and dimensionless spin 
parameters, the semi-major and semi-minor axis values 
(a, b) of the error ellipse on the sky, and the luminos- 
ity distance and angular resolution; the comparisons are 
shown in tables Hl lIIII 

The effect of choosing arbitrary coalescence times 
is illustrated in Fig. 2] where the distributions of lower 
bounds on the graviton Compton wavelength Xg are 
shown for fixed and random values of tc . It is clear 
from the graph that randomizing tc, leads to somewhat 
smaller lower bounds, with a tail at low values of the 
bounds, depicting the effect of signal loss in some of the 
cases. Fig. |4] also shows that using two LISA arm com- 
binations generally leads to improved bounds. 

In Fig. O we plot the distribution of the correlation 
coefficients between the massive graviton parameter /3g, 
and the two dimensionless spin parameters Xi (top panel) 
and the two masses rrii (bottom panel) for the 10^ -I- 
10^ Mq black hole case. The correlations are quite mild, 
with most of the values ranging between and 0.8, in 
contrast to the nonprecessing case where correlation 
coefficients larger than 0.9 were routine. This illustrates 
the strong decorrelating effect of the precessions. 

V. CONCLUSIONS 

In this paper we have studied bounds that can be 
placed on the mass of a hypothetical graviton using GW 
observations from the planned LISA mission, including 
spin precession effects. One possible extension of this 



work would be to include the effect of higher amplitude 
harmonics of the GW signal; in the non-spinning case, 
this is known to improve the accuracy of estimating pa- 
rameters, including distance and sky location (see, eg. 
(25I [26I \2^). and the graviton mass [6]. A final note: 
The inclusion of spin precessions has a significant com- 
putational cost in the parameter estimation procedure. 
Recently Kocsis et al. [2^ have developed a very efh- 
cient way, the harmonic mode decomposition, of decou- 
pling the several parameters that the signal depend on, 
according to their frequency "signature". This way the 
integrations for computing the elements of the Fisher ma- 
trices can all be done at once, lowering significantly the 
computational cost. It would be interesting to try to 
implement this decomposition on our code in the future. 
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mi {Mq) 


7712 (A^©) 


Ami /mi 


Am2 /m2 


Axi 


Ax2 


2a (arcmin) 


2b (arcmin) 


AQs (deg^) 


ADl/Dl 


3 X lO"" 


lO'' 


0.000667 


0.000541 


0.00157 


0.00306 


16.9 


7.3 


0.0233 


0.00240 






0.000387 


0.000314 


0.00130 


0.00176 


13.9 


8.4 


0.0245 


0.00236 


3 X lO'' 


lO'' 


0.00238 


0.00192 


0.00380 


0.00674 


32.3 


14.7 


0.0839 


0.00419 






0.00458 


0.00371 


0.00357 


0.00613 


23.8 


14.6 


0.0730 


0.00193 



TABLE L Comparison of median errors in selected parameters for two cases of 10* precessing binaries at z — 1. Semi- major 
axes of error ellipse on the sky parametrized by a and b; angular resolution is Af2s. Lang-Hughes results are quoted in the first 
line; our results are in italics in the second line. 



mi (Mq) 


m2 (Mq) 


Ami/mi 


Am^jm^ 


Axi 


Ax2 


2a (arcmin) 


26 (arcmin) 


AQs (deg") 


ADl/Dl 


3 X 10" 


10" 


0.00363 


0.00294 


0.00879 


0.0171 


92.5 


32.5 


0.656 


0.0126 






0.00225 


0.00182 


0.00671 


0.0140 


83.5 


49.0 


0.885 


0.0058 


3 X 10** 


lO'' 


0.0181 


0.0148 


0.0223 


0.0386 


142 


64.6 


1.65 


0.0193 






0.0129 


0.0103 


0.0130 


0.0290 


96.8 


58.4 


1.21 


0.0161 








TABLE H: Same 


as Table HI but for z 


= 3 






mi (Mq) 


m2 (Mo) 


Ami /mi 


Am2/m2 


Axi 


AX2 


2a (arcmin) 


2b (arcmin) 


AQs (deg") 


ADl/Dl 


3 X lO'' 


10" 


0.00811 


0.00658 


0.0193 


0.0359 


217 


95.8 


3.73 


0.0284 






0.00476 


0.00410 


0.0108 


0.0150 


201 


123 


5.22 


0.0145 


3 X lO'' 


10** 


0.0576 


0.0475 


0.0606 


0.107 


304 


139 


7.52 


0.0436 






0.0656 


0.0536 


0.0468 


0.112 


190 


116 


4.68 


0.0164 



TABLE HI: Same as Table H] but for z = 5 



